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Abstract 

A suitable correction of the Maxwell model brings to an enlarge- 
ment of the space of solutions, allowing for the existence of solitons 
in vacuum. We review the basic achievements of the theory and dis- 
cuss some approximation results based on an explicit finite-difference 
technique. The experiments in two dimensions simulate travelling 
solitary electromagnetic waves, and show their interaction with con- 
ductive walls. In particular, the classical dispersion, exhibited by 
the passage of a photon through a small aperture, is examined. 
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1 Introduction 

It is known, although rarely stated in books, that Maxwell equations in 
empty space do not admit finite-energy solitary waves among their solu- 
tions. By this we mean smooth electromagnetic waves with compact sup- 
port, travelling along straight-lines at the speed of light. One of the reasons 
to explain this fact can be attributed to the linearity of the Maxwell model, 
not allowing the due focussing of the signal on a constrained path. Such 
an impediment has stimulated the research of alternative nonlinear mod- 
els, mostly based on modification of the Lagrangian (see for instance [5]), 
granting the existence of soliton-like solutions and justifying in this way 
electromagnetic phenomena such as photons. 
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A step ahead in the comprehension of electromagnetic soHtary waves 
has been done in [11] and [12], where, in the framework of a self-consistent 
theory, explicit analytic solutions are carried out. The main argument is 
that the classical equations of the electromagnetism in vacuum are not 
capable to follow the evolution of finite-energy wave-fronts in the proper 
way (i.e., the one described for instance by the Huygens principle); this is 
due to the difficulty to impose both divE = and divB = at the same 
time on each point of the same front. As a matter of fact, if we define the 
wave- front as the surface envelope of the vector fields, one can easily show 
that requiring that the front evolves along its normal direction is incom- 
patible with the two free-divergence conditions taken together. Therefore, 
the relation divE = has been dropped, without this necessarily imply- 
ing the existence of point- wise electric charges. Here, we circumscribe our 
exposition by outlining some of the peculiarity of the new approach, while 
we refer to [12] for a detailed explanation of the origin and the possible 
physical implications. 

In particular, in this paper we devote our attention to some numerical 
simulations. It is soon evident that the removal of relation divE = is 
important for the construction of numerical algorithms in general, because 
most of the difhculties encountered in simulations are indeed consequence 
of the imposition of such a constraint. These include the efforts made to 
build approximation spaces satisfying some divergence-free conditions (see 
for instance [5], [1], [7], [S], [2], [I?]) or divergence corrections techniques 
(see for instance [TS], [TH], [12]). Some modifications of the Maxwell model, 
which may be in some way assimilated to the ones considered here, have 
been proposed at numerical level, in order to set up stable schemes (see [22] , 
p4 ), to handle boundary conditions (see [1], [3]), or for the treatment of 
wave propagation in linear non-dispersive lossy materials. These techniques 
are mainly adopted to overcome numerical troubles and are not intended 
to modifying the Maxwell model itself, as we are doing here. A survey of 
methods is given in [20], where the reader can also find an updated list of 
references. 

For simplicity, we will use explicit finite-differences (in particular the 
Lax-Wendroff method). The aim is to validate the theory in [12] with 
a series of simple experiments. To this purpose we follow the evolution 
of solitary waves and see how they react when encountering matter. We 
will mainly concentrate to the physical appearance and significance of the 
solutions. Thus, no theoretical issues will be discussed. Moreover, the 
discretization parameters will be small enough to get sufficiently accurate 
approximations, finalized to recover qualitative information. In no way we 
shall claim that our numerical approach is competitive, being conscious 
that the one proposed here can be certainly improved upon, as far as per- 
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formances versus costs are concerned. 

The author would hke to thank A. UgoUni (see [53]) for setting up the 
computatfonal code. 



2 The model equations 

By denoting with c the speed of hght, the classical Maxwell equations, in 
void three-dimensional space, in absence of electric charges, take the form: 

— = c^curlB (2.1) 

dt ^ ' 

divE = (2.2) 

— = -curlE (2.3) 

divB = (2.4) 

where the two fields E and cB have the same dimensions. 

In [12], part of the analysis is devoted to electromagnetic free-waves. 
These are solutions of the following set of model equations: 

(9E 

— = c^curlB - pV (2.5) 
r)R 

— = -curlE (2.6) 

divB = (2.7) 

E + V X B = (2.8) 

where p — divE, and V is a velocity vector field satisfying |V| = c. The 
field V is oriented as the vector field E x B. Note that relation (|2.8p 
is certainly satisfied for all electromagnetic waves where E is orthogonal 
to B and |E| = |cB|. These requirements are standard. In addition, we 
expect the new set of equations to admit a large space of solutions. We 
will check later that, contrary to the Maxwell model, solitary waves with 
compact support, as well as perfect spherical waves, are now among the 
solutions. On the other hand, if divE is relatively small (as it actually 
happens in many practical circumstances), then (|2.5p is as accurate as 
the corresponding standard Maxwell equation. Therefore, we expect the 
new model to be consistent with the existing ones, for a broad range of 
applications. 

It should be clear that the added term pV = cpJ, where |J| — 1, has the 
meaning of a current density fiowing in the direction of the rays at speed c 
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(Ampere term). We assume this to be true even if there is no presence of 
point- wise electric charges. As documented in [12], such a current is part of 
the wave itself. Moreover, equation (|2.8p actually characterizes free-waves, 
since its says that there are no external 'forces' acting on the wave-fronts, 
as a result of external perturbations. It is a sort of Lorentz law (where 
the force turns out to be zero), in which moving charges are replaced by a 
balance of pure vector fields. We will see later, in section 4, how equation 
(|2.8p should be modified in order to handle situations in which the wave 
interacts with matter. 

It is important to observe that a conservation principle holds. In fact, 
by taking the divergence of equation (j2.5p , we come out with the following 
continuity equation: 

^ = -cdiv(pJ) (2.9) 

In addition, using that V is orthogonal to both E and B, one can easily 
prove the classical Poynting relation: 

^^(|E|^ + |cB|2) = -c^div(ExB) (2.10) 

Such properties confirm that what we are doing has physical meaning. 
It can be proven that the electromagnetic free-waves described by the new 
set of equations are perfectly compatible with the Huygens principle and 
the eikonal equation (see for instance [Q). In fact, if, for example, V is 
irrotational (i.e.: the rays are straight-lines), then V = V^" for some poten- 
tial 4", so that the relation |V| = c is equivalent to IV^*! = c, which is the 
eikonal equation. Moreover, ()2.5p follows from minimizing the standard 
Lagrangian of classical electromagnetism, after imposing the constraint 
A = $J to the electromagnetic potentials A and $ (see [13] and [12], 
section 2.4). The covariant version of the equations and their invariance 
under Lorentz transformations are also discussed in [12], chapter 4. 



3 The discretization method for 2-D soUtons 

Explicit expressions of wave-packets, moving at speed c and solving the 
system of equations (|2.5p - p.8|) . are available. The electromagnetic signals 
carried by these fronts can be very general. For the sake of simplicity, only 
examples in two dimensions will be examined here. This does not obscure, 
however, the general validity of the formulation. 

We start by discussing a simple case. In the 3-D space with cartesian 
coordinates {x,y, z), we can define the following fields: 

E = (^cfix)gict-z), 0, O) 
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B = (O, fix)g{ct-z), O) V = (0, 0, c) (3.1) 

Now, the continuity equation (|2.9p becomes the simple transport equa- 
tion: 

In practice, the example is two-dimensional since there is no dependency 
on the variable y. The fields describe the evolution of a solitary wave- 
packet shifting along the axis z. This is a group of parallel transversal 
fronts, travelling in the same direction V and modulated in intensity by 
the function g. It is important to remark that the fields E and B satisfy 
Maxwell equations (including p.2|) ) only if the function / is constant (this 
check can be done by directly differentiating the expression in (|3.1|) and 
substituting in the model equations). This corresponds to a plane wave of 
infinite energy, which is the only Maxwellian front allowed by the above 
setting. If we want the energy to be finite, the only possible choice is / = 0, 
implying E = 0. The situation is different if we take into account equation 
(|2.5p . Now, as the reader can easily check by direct computation, the fields 
in Ij3.ip are solutions, and the two functions / and g can be completely 
arbitrary. In particular, / and g may have compact support. Note that, 
for simplicity, the wave is unbounded in the y direction, hence, in this case 
the energy of the wave-packet is not finite. Examples where the fronts are 
also bounded in the y direction can be easily constructed (recall to enforce 
the condition divB = 0), but the problem becomes three-dimensional and 
here we would prefer to avoid this complication. 

Let us discuss a specific example. We start by fixing a point (xq,zq). 
At time t — 0, we can assign the distribution: 



(0^(2:)(/)^(z), 0, O) with 



4>t{x) = 1 + COS ^^^^ —TTJ, 4>aiz) = 1+cos^^^ — —Trj (3.3) 

inside the rectangle centered at {xq,zo) with sides of length 2t and 2a. 
Outside the same rectangle E is required to vanish. Accordingly, field cB 
is supposed to be orthogonal to the plane {x, z) and with the same intensity 
as E. When time passes, this signal-packet moves at speed c in the direction 
of the axis z: 

E = (cj)r{x)(l>„{ct-z), 0, O) (3.4) 

We do not need to use any approximation technique in order to follow 
the evolution of such a soliton, since the corresponding analytic expression 
exactly solves the set of equations (j2.5|) - (|2.8|) . However, in view of more 
complex applications, let us introduce the Lax-Wendroff method. More 
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precisely, we work with its 2-D version applied to a system of pure hyper- 
bolic equations, as described in [21], section 4.9. 

Let us suppose that, for some real constant dx d matrices and M^, 
our system can be written in the following form: 



du du du 

— = M,— + M,— 
ot ax oz 



(3.5) 



where the d-components vector u depends on x, z and t. 

We denote by {xt, Zj) the nodes of a uniform grid (of step-size h) on the 
plane {x,z). The time-discretization parameter is denoted by At. Based 
on the stencil shown in figure 1, we pass from the approximated solution 
u'^ at time tk to the one at time tk+i ~ tk + At. 
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A = {xi^i, 

B = ixi-i,Zj) 
C = (xi-i, Zj+i) 

D = {Xi, 2j_l) 

E = {xt, Zj) 
F = {x^,Zj+i) 
G = (xi+i, 
H = {xi+i,Zj) 

I = {x,+i,Zj+i) 



Figure 1: Stencil for the Lax-Wendroff method. 



This is done according to the scheme: 

u^+^{E) = \l -\^{Ml + Ml)\u''{E) + \\A'L^{I + \M^)u^{H) 

- iAM,(/-AM,)u'=(B) + \\M,{I+\M,)m''{F) - \\M,{I-\M,)vi\D) 
+ \\^{M,M, + M,M,)(u'=(/) - u'=(G) - m\C) + u^(A)) (3.6) 

where A = l^t/h and / is the identity matrix. 

Let us examine the case of a wave with the electric field laying on the 
plane (x, z). Therefore, the magnetic field, which is orthogonal to that 
plane, has only one component different from zero. The condition (|2.4|) is 
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always satisfied, since our functions do not depend on the variable y. We 
organize the vector u in the following way (d = 3): 



u ^ {El, cB2, E3) 

For the classical Maxwell system (see (|2.ip and (|2.3p ') it is enough to define 
the matrices: 





= c = -c (3.7) 



in order to recover (j3.5p . It is important to remark, however, that the condi- 
tion (|2.2p is not enforced (this should be done at the initial time, although, 
as we mentioned, fronts having divE = and divB = are extremely rare 
and out of interest for our applications). Note that the eigenvalues of both 
matrices are ±c, confirming that the system is of hyperbolic type. 

Regarding the modified system (see ()2.5p and (12. 6p ). the above matrices 
have to be replaced by: 





iV4 = c = -c (3.8) 



where V = (Vi,0, V3) is a velocity vector orthogonal to E = (£'i,0, £^3), 
such that Vi = c^. Moreover, we require relation (j2.8p to be verified. 

In other terms, V must be defined as: c(E x B)/|E x B|. Note that this 
time we do not need to satisfy (|2.2p . Finally, let us observe that, for both 
the matrices in (|3.8p . the modulus of the maximum eigenvalue turns out 
to be equal to c. 

The matrices in p.8p are not constant. Nevertheless, we will apply the 
Lax-Wendroff method to the new system also in this circumstance. There- 
fore, when implementing (|3.6p at the fc-th step, the entries of the matrices 
Mx and Mz will be evaluated at the central point of the stencil of figure 
1. We are aware of the fact that this correction may deteriorate a bit the 
convergence properties of the method. 

In order to guarantee stability, we must impose the following constraint 
(CFL condition): 

In our experiments we always used the maximum At allowed by p.9p . i.e.: 
At = h/2cV2. With this choice we never encountered problems concerning 
stability. 
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We start by testing the discretization method on the sohton correspond- 
ing to the initial conditions in Ij3.3p . with a = t. The initial condition for 
the discrete scheme is obtained by interpolation at the grid-points. It is 
interesting to note that, with this setting, the second component of E re- 
mains zero during the evolution. This is true both for the exact solution 
(see ^6A\i ) and the discrete one. As a matter of fact, the scheme (I3.6P 
follows from repeated applications of the three matrices /, and Mz ■ In 
this situation we have Vi = and = c. Every time we apply the matrix 
Mx to a vector of the form *, 0), we get a vector of the same form, since 
El = 082- The same is true regarding the matrix M^. It is also interesting 
to note that relation (|2.8p is preserved for any k. 



Figure 2: Level sets of the scalar quantity jE|, at different times, for a free 
wave-packet shifting parallel to the z axis (from right to left). The vector E 
is orthogonal to the direction of motion individuated by V. The vector B is 
orthogonal to the page. 

We can see in figure 2 the moving soliton at different times. The dis- 
cretization parameters have been taken small enough to get rid of numerical 
disturbances. In other words, the parameter a is large enough to allow the 
inclusion of a sufficiently high number of grid-points in the support of the 
wave (there arc 13 x 13 nodes in this example). It is interesting to remark 
that the scalar function Ei satisfies the wave equation: 






dEi _ ^dEi 



(3.10) 
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Note instead that the wave equation — c AE does not hold, since 

in p.lOp we are omitting the derivatives with respect to the variable x. On 
the other hand, we are not approximating a trivial transport problem for 
the unknown E\, but a coupled hyperbolic system involving both E and 
B. This is not a trivial remark, but a rather important achievement. 




Figure 3: Level sets of jE|, at different times, obtained by implementing the 
scheme (|3.6p with the matrices given in (|3.7p . 

Some longitudinal numerical dissipation is visible. It can be suitable 
reduced by using more appropriate numerical techniques. It is important 
however to remark that no numerical viscosity is introduced in the direction 
transverse to the motion, so that the solitary wave proceeds maintaining a 
well defined width. This is quite an important result. Suitable Sommerfeld 
type boundary conditions are assumed at the outflow boundary, in order 
to allow a smooth exit of the wave from the computational domain. 
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A completely different situation (see figure 3) is observed when using 
the same initial datum (|3.3p and evolving the wave with the help of the 
two linear Maxwell equations (|2.ip and (|2.3p . using the matrices in (|3.7I) . 
In this case the condition (|2.2p is not verified, since such a constraint is not 
satisfied by the initial datum and it is not enforced by the numerical scheme. 
As the reader can see, part of the solution is diffused all around, despite 
the fact that, in the initial dislocation, the vectors V are all oriented in the 
horizontal z direction. Surprisingly, the most of the energy does not move at 
all. Our guess is that the condition divE ^ generates a kind of stationary 
charge. Of course, this behavior is not physically acceptable. Nevertheless, 
there is no way to modify the initial datum, maintaining / and g with 
compact support. In fact, as mentioned before, one can verify that in (|3.1[) 
the only possible initial condition compatible with (|2.2[) is E = 0. As far 
as solitons are concerned, this indicates that the poor performances of the 
classical equations of electromagnetism are mainly due to the total lack of 
initial compatible conditions. 




Figure 4: The first picture shows the level sets of |E|, for some initial datum 
of the form given by (|3.11|l . The two successive enlargements show the initial 
distributions of E and V in the plane (a;, z). 



A less trivial example is obtained by taking an initial datum of the 
form: 

E = (e^{x,z), 0, E^{x,z)^ (3.11) 

whose shape is given in figure 4. The magnetic field is along the direction 
of the y-axis and such that |E| = |cB|. The corresponding velocity field V 
is radially distributed with respect to a prescribed point. 
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Figure 5: Evolution of the wave corresponding to the initial condition of figure 
4, as obtained by running the numerical scheme (|3.6p with a relatively small 
discretization parameter. 

With this setting, during the evolution one gets a semi-cylindrical wave, 
simulated by the numerical scheme as shown in figure 5. There is a decay 
in intensity of the electromagnetic fields, due to the fact that the signal 
is spread out on front surfaces having increasing magnitude. This is in 
agreement with the equation of continuity and the energy preservation rules 
(see (|2.9p and (|2.10p ). automatically imposed by the governing equations. 
The approximation is not excellent, since it is polluted by a secondary wave 
of small amplitude evolving backwards. We think that this is due to the 
poor performances of the Lax-Wendroff scheme. In fact, in the non-convex 
part of the support of the wave some oscillations are generated. They 
modify the sign of the electric field, but not the one of the magnetic field, 
so that the associated velocity field V (proportional to the Poynting vector 
E X B) locally changes orientation. This anomalous part is then correctly 
simulated by the numerical method, but, unfortunately it develops along 
a wrong path. This is a minor effect that should be however kept in mind 
when designing an alternative code. 

After testing our algorithm on these simple examples, we are now ready 
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to deal with a more serious situation. 



4 Constrained waves 

A generalization of the model is obtained by modifying equation (|2.8p in 
the following way (see [T^], chapter 5): 



where the constant > is a charge divided by a mass, p — divE and p is 
a kind of pressure. As customary, the substantial derivative G — -^V is 
defined as -^V = ^V+(V-V)V, so that G turns out to be an acceleration. 
Basically, if V is the normalized vector field tangent to the light rays, then 
the field G gives a measure of their curvature. 

Equation ()4.1|) . even if there are no classical moving charges, can be 
assimilated to the Lorentz law. If G and p are zero, then one is dealing 
with a free electromagnetic wave, and the corresponding rays are straight- 
lines. Actually, equation (|2.8[) . where G = and p = 0, says that the 
evolution of the wave is free from constraints, that is, there are no external 
factors ('forces') acting on it. When, for some reasons, (|2.8p is not satisfied, 
then pressure develops. At the same time -^V is different from zero, so 
that V changes direction (the rays are curving) and the electromagnetic 
wave-fronts locally follow the evolution of the new normalized Poynting 
vector J = V/c. Thus, when G is different from zero, the wave is no longer 
free, and, following [12], it will be called constrained wave. Such a situation 
happens when the wave is subjected to external electromagnetic fields. This 
is true for instance during the interaction with matter at atomic level, as 
in reflection or diffraction. Therefore, our model equations are now able 
to provide the coupling between the curvature of the rays and the motion 
of the wave-fronts. Note that, during the change of trajectory of the rays, 
depending on the context, the polarization may also vary. 

It is evident that equation ()4.ip is inspired by the Euler equation for 
inviscid fluids (with an electromagnetic type forcing term). Put in other 
words, the vector field V evolves according to the laws of fluid dynamics. 
At the same time, it carries, on the transversal manifold, the two orthogonal 
fields E and B. In this peculiar way, the system of equations (|2.5p - (|2.6p - 
(P?7)l - ((i7T|) describes the evolution of the triplet (E, B, V). We remark that 
our model is a system of vector equations. Therefore, although our anal- 
ysis is based on simple examples, it has nothing to do with the equations 
modelling scalar solitons (see for instance [10] )• 

The exact physical meaning of (|4.ip is carefully explained in [12] , where 
an equation of state involving p, that uses the curvature of the space-time 




(4.1) 
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geometry, is also introduced. For the sake of simplicity, we do not discuss 
further the reasons of this construction and we refer to [12] for the details. 
Numerical experiments on constrained waves, concerning photons trapped 
in toroid shaped regions, are illustrated in [5]. 

5 Some experiments on diffraction 

Our next step is to analyze the interaction of an incoming soliton with a 
perfectly conducting plane wall. We would like to discuss some experiments 
concerning diffraction, so that there is a hole in the wall and our soliton will 
be forced to pass through it. The barrier is simulated by imposing suitable 
boundary conditions, so that the impact is immediate and discontinuities 
are generated. 




Free wuve 



Figure 6: An example of diffraction. A soliton is scattered when passing through 
an aperture smaller than its size, producing a diffusive effect. 



According to figure 6, we expect the soliton to behave as a free- wave 
(hence G = and p = in (|4.ip ') before reaching the aperture. The 
reaction of the wall brings to an instantaneous diffusion of the signal, as 
also observed in real-life experiments with fluids or electromagnetic waves. 
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After the passage, the sohton continues his journey behaving as a free wave, 
following however non-parallel characteristic straight-lines. 

Assuming that the path of the incoming wave is orthogonal to the 
boundary, there we impose V = (no slip condition). Hence, for an 
instant, the last term on the right-hand side of equation (|2.5p disappears. 
This way of imposing boundary conditions is suggested in |12| . chapter 5, 
where the following Bernoulli type equation is deduced: 

|^(p+fl|Vf) ^ -fllVfdivV (5.1) 

Heuristically, when V goes to zero at the boundary, the pressure raises 
from zero to infinity, in order to maintain the energy balance. One can 
actually check that divV > after the impact, bringing to a diffusive 
effect (with no dissipation of energy, however). As we shall see from the 
experiments, these assumptions seem to be correctly posed. 




Figure 7: Evolution of a solitary wave-packet passing through the hole of a 
conductive wall slightly smaller than its width. 



From the numerical point of view the problem is stiff, therefore some un- 
avoidable perturbations will be observed. A more realistic behavior could 
be obtained by smoothing the effects of the boundary. This could be sim- 
ulated by slowing down the incoming soliton by varying V in a continuous 
way. We did not investigate however this possibility. 
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Figure 8: Evolution of a solitary wave-packet passing through the hole of a 
conductive wall, having smaller size if compared to the one considered in figure 
7. 



In figures 7 and 8, we can see the evolution of the norm of the electric 
field. The width of the incoming soliton includes 32 grid points, while the 
size of the hole corresponds to 28 grid points in figure 7 and 26 grid points 
in figure 8. As expected, the no slip condition at the boundary produces 
a diffusive behavior. After hitting the obstacle, the development of the 
soliton is similar to that of figure 5 (with much less emphasis, the fronts 
assume a circular asset, and we know from the experiments of section 3 
that the scheme can support this situation). There is no loss of energy 
or introduction of transversal dissipation. The spreading is a consequence 
of the new geometrical distribution of the rays, following non-parallel tra- 
jectories (in other words, as in figure 6, the stream-lines associated to the 
velocity vector V are diverging). 

There are little wiggles scattered all around. These are mainly due to 
the exceeding part of the wave, reflected back by the wall. Therefore, these 
disturbances might be (in part) physically consistent. Meanwhile, they can 
also point out the limit of the algorithm and suggest better ways to handle 
boundary conditions (especially for the portion involved in the reflection 
process). 
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Finally, in figure 9, we tested the case when the incoming soliton is not 
symmetrically centered with respect to the aperture. Of course, we could 
show a multitude of other scattering experiments, regarding wave-packets 
of arbitrary shape interacting with a variety of obstacles. However, the 
purpose here was to provide a general qualitative insight. 



Figure 9: Evolution of a solitary wave-packet passing through the hole of a 
conductive wall. The aperture is now asymmetrically placed with respect to 
the trajectory of the wave-packet. Some minor disturbances are present. To 
some extent, they can be physically acceptable, but they are also due to the 
diffusive effects of the Lax-Wendroff method, contrasted by the sharp boundary 
constraints. 



These experiments can certainly be ameliorated by using more appro- 
priate approximations methods. Rather than accurate results, the interest 
here was to show that, with the help of the model equations (|2.5|) - (|2.6p - 
(|2.7p - (l4.ip . an electromagnetic wave can be treated exactly as a material 
fluid. We would like to remark once again that, although the pictures pre- 
sented here are related to the evolution of the scalar quantity |E|, we are 
not solving a trivial scalar equation. Our solutions are electromagnetic lo- 
calized emissions, presenting both the features of waves and particles (pho- 
tons). Their internal structure is perfectly determined in terms of electric 
and magnetic fields. In the more complex and realistic 3-D case, they can 
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even change polarization, reacting in this way like true electromagnetic en- 
tities, but with energies concentrated in finite regions of space. We believe 
that these peculiarities are important both for physical implications and 
technical applications. 
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